clear
set more off
set type double
graph drop _all
program drop _all
set scheme s2mono
//log using "\\C:\Users\XUE XINDONG\Desktop\meta_educ_health_part1.smcl", replace
use "C:\Users\XUE XINDONG\Desktop\educ_health\meta_educ_health(20190928).dta",

******************************************************************************
*DEFINITION OF BASIC VARIABLES
******************************************************************************
gen tstat=tstatistics
gen negative = 1
replace negative = -1 if BadHealth == 1 & BadEduc == 0
replace negative = -1 if BadHealth == 0 & BadEduc == 1
replace tstat = tstat*negative
gen df = n-k1
gen r = tstat/sqrt(tstat^2+df)
gen varR = (1-r^2)/df
gen seR = sqrt(varR)
gen pcc = r
gen abststat = abs(tstat)
gen abspcc = abs(pcc)

******************************************************************************
* Distribution of t-stats and PCCs before kicking out outliers
******************************************************************************
summ tstat, detail
summ pcc, detail
summ df, detail 


******************************************************************************
*Kicking out outliers
******************************************************************************
reg pcc seR
predict resid, residual
egen stdr=std(resid)
gen absr=abs(stdr)
drop if absr>3.5
drop resid stdr absr


******************************************************************************
* Characteristics of t-stats and PCCs of final data set
******************************************************************************
summ tstat, detail
histogram tstat,bin(80)

histogram tstat, bin(25) name(histTSTA1)
gen tstat1 = (tstat<-2)
gen tstat2 = (tstat>=-2 & tstat<=2)
gen tstat3 = (tstat>2)

ge tstatn=0
replace tstatn=1 if tstat1==1
replace tstatn=2 if tstat2==1
replace tstatn=3  if tstat3==1
ta tstatn

// NOTE: A lot of insignificant t-stats
sum tstat1 tstat2 tstat3



summ pcc, detail
histogram pcc, bin(80) name(histPCC)

histogram pcc if PhysicalHealth1==1, bin(80) name(histPHPCC)
histogram pcc if Mentalhealth1==1, bin(80) name(histMHPCC)
histogram pcc if Generalhealth1==1, bin(80) name(histGHPCC)

sum PhysicalHealth1 Mentalhealth1 Generalhealth1

******************************************************************************
* Relationship between t-stats and PCCs
******************************************************************************
// This gives the correlation between pcc and tstat. It's not real high (0.61).
// Next I graph the two to show the positive correlation.
corr pcc tstat
graph twoway (scatter pcc tstat, msize(*0.5) msymbol(Oh) name(PCCtstat)) (lfit pcc tstat) 
summ df, detail

******************************************************************************
* Relationship of PCCs over time
******************************************************************************

bysort ID: egen PCCmed = median(pcc)
graph twoway (scatter PCCmed PubYear, msize(*0.5) msymbol(Oh) name(PCCtime21)) (lfit PCCmed PubYear) 


******************************************************************************
**CALCULATING THE MEAN EFFECT WITHOUT CORRECTING FOR PUBLICATION BIAS
******************************************************************************


******************************************************************************
* Calculating study weights based on number of estimates per study
******************************************************************************
bysort ID: egen numberests = count(ID)
gen weight = 1/numberests


******************************************************************************
* FIXED EFFECTS
******************************************************************************

gen feprecisionR = 1/seR
gen fetstatR = r/seR


******************************************************************************
* Equal weight to each estimate (Column 1)
******************************************************************************
regress fetstatR feprecisionR, noc vce(cluster ID)


******************************************************************************
* Equal weight to each study (Column (2)
******************************************************************************
regress fetstatR feprecisionR [pweight = weight],  noc vce(cluster ID)


******************************************************************************
* RANDOM EFFECTS 
******************************************************************************
//metan r seR, random
metareg r seR , wsse(seR)
scalar tau2 =  e(tau2)
gen revarR = varR + tau2
gen reseR = sqrt(revarR)


******************************************************************************
* Correcting for heteroskedasticity
******************************************************************************
gen reprecisionR = 1/reseR
gen retstatR = r/reseR
gen repubbiasR = seR/reseR


******************************************************************************
* Equal weight to each estimate (Column 3)
******************************************************************************
regress retstatR reprecisionR ,  noc vce(cluster ID)


******************************************************************************
* Equal weight to each study (Column 4)
******************************************************************************
regress retstatR  reprecisionR [pweight = weight] , noc vce(cluster ID)



******************************************************************************
* DETECTING PUBLICATION BIAS
******************************************************************************

******************************************************************************
* Funnel Plot for individual estimates
******************************************************************************
// FUNNEL GRAPH
metafunnel r seR, ylabel(#8) name(funnelA) //name(funnel, replace) nolines
metafunnel r seR if seR < 0.4, ylabel(#8) name(funnelAA) //name(funnel, replace) nolines
//graph save funnel1, replace
//NOTE: No evidence of asymmetry in estimates, hence no evidence of publication bias
//NOTE: But much evidence of random effects!

******************************************************************************
* Funnel Plot for study means
******************************************************************************
//FUNNEL GRAPH
by ID, sort: egen meanr = mean(r)
by ID, sort: egen meanseR = mean(seR)
metafunnel meanr meanseR, ylabel(#8) name(funnelB)  //name(funnel, replace) nolines


******************************************************************************
* FAT/PET
******************************************************************************

******************************************************************************
* FIXED EFFECTS
******************************************************************************

******************************************************************************
* Equal weight to each estimate (Column 1)
******************************************************************************
regress fetstatR feprecisionR,  vce(cluster ID)


******************************************************************************
* Equal weight to each study (Column (2)
******************************************************************************
regress fetstatR feprecisionR  [pweight = weight],  vce(cluster ID)


******************************************************************************
* RANDOM EFFECTS 
******************************************************************************

******************************************************************************
* Equal weight to each estimate (Column 3)
******************************************************************************
regress retstatR  reprecisionR repubbiasR,  noc vce(cluster ID)


******************************************************************************
* Equal weight to each study (Column 4)
******************************************************************************
regress retstatR  reprecisionR   repubbiasR [pweight = weight] , noc vce(cluster ID)



******************************************************************************
* *PEESE
******************************************************************************

******************************************************************************
* FIXED EFFECTS
******************************************************************************

******************************************************************************
* Equal weight to each estimate (Column 1)
******************************************************************************
regress fetstatR feprecisionR seR,  noc vce(cluster ID)


******************************************************************************
* Equal weight to each study (Column (2)
******************************************************************************
regress fetstatR feprecisionR seR [pweight = weight], noc vce(cluster ID)


******************************************************************************
* RANDOM EFFECTS 
******************************************************************************
ge  repubbiasRR = varR/reseR

******************************************************************************
* Equal weight to each estimate (Column 3)
******************************************************************************
regress retstatR  reprecisionR   repubbiasRR,  noc vce(cluster ID)


******************************************************************************
* Equal weight to each study (Column 4)
******************************************************************************
regress retstatR  reprecisionR repubbiasRR  [pweight = weight] , noc vce(cluster ID)


//log close
 